#include <stdio.h>
#include <math.h>
//注意：本次拟合因为精度要求，使用了long double类型的函数
//然而win环境下DevC++并不支持Lf的输出，所以做了一些针对性处理 
//最终代码只能用gcc编译 

int main() {
    long double x[5] = {2, 3, 4, 5, 6};
    long double y[5] = {10, 8, 1, 8, 10};
    long double N = 5.0; // 样本个数，为了避免影响后续乘除选用long double形式
    long double a = 5; // 初始化a
    long double b = 0; // 初始化b
    long double c = 0; // 初始化c
    // 不断调整迭代a、b、c的大小，直至函数走到最低点，解出y=ax^2+bx+c
    long double alpha = 0.001; // 初始化学习率，即梯度下降的步长
    //__mingw_printf("%.12Lf",x[0]); 在重复输出0后尝试寻找原因，发现是编译器的问题 

    // 梯度下降部分
    int max_iterations = 1000; // 最大迭代次数，达到最大迭代次数时停止迭代
    for (int iter = 0; iter < max_iterations; iter++) {
        long double dJ_da = 0; // 对a求偏导，初始化
        long double dJ_db = 0; // 对b求偏导，初始化
        long double dJ_dc = 0; // 对c求偏导，初始化
        
        for (int i = 0; i < (int)N; i++) {
            dJ_da -= 2.0 * x[i] * x[i] * (y[i] - (a * x[i] * x[i] + b * x[i] + c)); // 对a求偏导(还没除以N)
            dJ_db -= 2.0 * x[i] * (y[i] - (a * x[i] * x[i] + b * x[i] + c)); // 对b求偏导(还没除以N)
            dJ_dc -= 2.0 * (y[i] - (a * x[i] * x[i] + b * x[i] + c)); // 对c求偏导
            //这里利用了分别求导可以累加的特点简化原方程的运算 
            __mingw_printf("%Lf %Lf %Lf\n",dJ_da,dJ_db,dJ_dc);
        }
        dJ_da /= N;
        dJ_db /= N; // 在循环外做最后一次处理，得到偏导，此处不能用int
        dJ_dc /= N;
        a -= alpha * dJ_da; // 步长乘以偏导为在a轴移动的距离
        b -= alpha * dJ_db; // 在b轴移动的距离
        c -= alpha * dJ_dc; // 在c轴移动的距离
    } // 循环，得到a，b，c

    __mingw_printf("The fitted curve is y = %Lf x^2 + %Lf x + %Lf\n", a, b, c);

    return 0;
}
//在解决了输出问题后发现迭代出的a、b、c过于大了，迭代处于发散的状态 
//所以调整了学习率，之前的学习率过大导致函数不能收敛 
//学习率可以动态调整